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We present an approach for flux analysis in process algebra models of biological systems. We per- 
ceive flux as the flow of resources in stochastic simulations. We resort to an established correspon- 
dence between event structures, a broadly recognised model of concurrency, and state transitions of 
process models, seen as Petri nets. We show that we can this way extract the causal resource depen- 
dencies in simulations between individual state transitions as partial orders of events. We propose 
transformations on the partial orders that provide means for further analysis, and introduce a software 
tool, which implements these ideas. By means of an example of a published model of the Rho GTP- 
binding proteins, we argue that this approach can provide the substitute for flux analysis techniques 
on ordinary differential equation models within the stochastic setting of process algebras. 



1 Introduction 



Computer science methodologies are now gaining increasing attention in systems biology, where they 
are used together with techniques from applied mathematics and physics to provide insights into the 
workings of biological systems. In such a setting, it is often desirable to model smaller biological systems 
as sets of differential equations. Being equipped with well understood analysis techniques, differential 
equations provide the deterministic approach (as opposed to the stochastic approach) for modelling and 
simulating biological systems. Differential equation models, however, are inherently difficult to change, 
extend and upgrade when modification to the model structure is required. In addition, within the last few 
years, a general consensus has emerged that noise and stochastic effects, which are not directly captured 
by differential equations, are essential attributes of biological systems (see, e.g., ||25]| ). 

'Tailored to fit' constructs of computer science languages provide the means for building models 
of larger systems, especially when it is required to refine, modify or compose the models, for example 
when new data is acquired. Various languages with stochastic simulation capabilities based on, for 
example, process algebras IIH |6l [HI |23l, term rewriting (see, e.g, S [131) and Petri nets (see, e.g., 
|[T6ll26]| ) have been proposed. In particular, process algebras, which were originally introduced to study 
the properties of complex, reactive, concurrent systems are well suited for modelling biological systems. 
By using these languages, we can build models to describe the dynamics of the biological systems, which 
are inherently complex and massively parallel, and run stochastic simulations. The time series data 
resulting from these stochastic simulations are then used to observe the emergent behaviour, together 
with other techniques that are borrowed from theoretical computer science to explore the model for 
hypothesis generation. However, these techniques are still underdeveloped in comparison to the rich 
arsenal of differential equation analysis techniques that have their roots in Newton's physics. There 
is now a pressing need for the adjustment of the computer science analysis techniques on models of 
biological systems and on stochastic simulations with them so that biologically meaningful observations 
can be made. 

Along these lines, we present an approach for the analysis of stochastic simulations with process 
algebra models. The stochastic semantics of process algebras are conveniently given by continuous time 
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Markov chains (CTMC). Each simulation is a trajectory of a CTMC, which is a sequence of computations 
of the underlying transition system, given by the process algebra model. The CTMC semantics, thus, 
imposes a total order on the transitions of a simulation trajectory that is emphasised by the unique time 
stamps of the individual transitions. In this respect, a simulation on a model can be seen as reduction 
of a complex structure, i.e., the model, into a simpler structure, i.e., the simulation trajectory. However, 
during this reduction some of the information on the model is lost, and others are made implicit. The idea 
here is to recover this implicit information: when these transitions are inspected from the point of view 
of their dependencies on one another, it is possible to relax the total order of the transitions into a partial 
order. We can then use this partial order as a representation of causal dependencies in the simulation and 
process it to observe the flux in the system with respect to the flow of the resources during simulation 
from one transition to another. 

Partial orders reflecting dependencies in computations have been studied as non-interleaving models 
of concurrency: when there are no resource conflicts between two partially ordered events, they can take 
place in parallel, and their common predecessors and successors provide a synchronisation mechanism. 
Event structures reflect this view [ 21i 3.1} : in an event structure dependencies of events of a concurrent 
system are given with a partial order. In ifTSI . for sequences of computations in Petri nets, we have given 
an algorithm for extracting partial orders that exhibit event structure semantic^ In the following, we 
carry these ideas to CTMCs of process algebra simulations by tracing the consumption and production 
effect of each event. This way, we extract the information on the flow of resources in simulations as 
partial orders. We then apply transformations to these partial orders for flux analysis. We illustrate these 
ideas first on a simple example, and then on a published model of Rho GTP-binding proteins |[T4l l5l. 

We restrict the presentation in this paper to the models that are given as sets of chemical reactions, 
and we use the stochastic n calculus and its implementation language SPiM |22|. However, we believe 
that the ideas presented in this paper are applicable for a more general class of models and to a broader 
spectrum of modelling languages. 

2 Processes of Biology 

In this section, we describe how stochastic n calculus can be used to model biological systems, e.g., as 
systems of chemical reactions |5| or as systems of entities with more complex structures, e.g., polymers 
m. In stochastic n calculus, the basic building blocks are processes which are defined as follows. 

Definition 1 /J'/ Syntax of the stochastic n calculus: processes range over P,Q,. . . Below fn{P) denotes 
the set of names that are free in P. 
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The algorithm in [18| is presented on multiplicative exponential linear logic encodings of Petri nets. 
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Expressions above are considered equivalent up to the least congruence relation given by the equiv- 
alence relation = defined as follows. 

P I () = P 

P I Q = Q I P 

P I (Q I R) = (P I Q) I R 

X(m) = P X(n) = P{m:=n] 

new X () = 

new X new y P = new y new x P 

x^fn(P) new x (P I Q) = PI new x Q 



While implementing models in the stochastic K calculus, we follow the abstraction of the biochem- 
ical species as entities that can change state. The state changes can be spontaneous or result from the 
interactions of the species with other species. We denote different states of the species with processes. 
Because the state of a species is a process with different state transition capabilities, a process P can 
choose stochastically between zero or more alternative behaviours. In the language of SPiM, we write 
the choice of processes as do PI or ... or PN. A choice of only one process is written as PI, 
while the empty choice is written as () . A parallel composition of N processes is written as PI I ... 
I PN. A process can also be given a name X with parameter m, written let X(m) = P. By using these 
constructs, we can compose processes to construct models and incrementally extend them to build larger 
models when desired. 

The reduction rules of the calculus are as follows. 

Definition 2 /5/ Reduction in the stochastic K calculus. 

(1) do delay@r; P or ... — ^ P 
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A species can spontaneously change state after a delay or as a result of a stochastic interaction 
with another species. The stochastic behaviour of the processes is delivered by the rate of the delay or 
interaction actions that the processes can perform: a process can perform a delay at rate r and then do P, 
written delayOr ; P. The rate r is a real number value denoting the rate of an exponential distribution. A 
process can also send a value n on channel x with weight wi and then do Pi , written ! x (n) *wl ; PI, or it 
can receive a value m on channel x with weight W2 and then do P2, written ?x (m) *w2 ; P2. With respect 
to the reduction semantics of SPiM, these complementary send and receive actions can synchronise on 
the common channel ;c and evolve to PI I P2{m := n}, where m is replaced by « in process P2. This 
allows messages to be exchanged from one process to another. Each channel name x is associated with 
an underlying rate given by p {x) . The resulting rate of the interaction is given by p {x) times the weights 
w\ and W2. If no weight is given then a default weight of 1 is used. The operator new xOr : t P creates a 
fresh channel x of rate r to be used in the process P, where t is the type of the channel. When a process is 
prefixed with the declaration of a fresh channel, that channel remains private to the process and does not 
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new bindOr 1 : chan ( chcoi) 

let A() = ( new e@r2 : chan new bindOrl : chanO 

!bind(e); Ab(e) ) let A() = !bind; AB() 

and Ab(e:chan) = !e; A() and AB() = delay(§r2; ( A() I B() ) 

let B() = ( ?bind(e); Bb(e) ) and B() = ( ?bind() ; () ) 

and Bb(e:chan) = ?e; B() 

Figure 1: Two different encodings in S Pi M of a model of the association of the species A and B to form 
a complex, and their dissociation: on the left, we have a model and its graphical representation, where 
the usage of the private name e for the bond between the species Ab and Bb permits the monitoring of 
the individual species in the complex after the association. The model on the right abstracts away from 
this information within a simpler representation. 

conflict with any other channel. By using these constructs, we can describe the dynamics of a biological 
system as a model in the language of SPiM. 

Example 3 Consider the biochemical species A and B which can associate with rate r\ to form a complex 
and dissociate with rate rj. In Figure^ we depict two different encodings of this model in SPiM. 

For an in depth exposure to the syntax and semantics of the stochastic K calculus that is implemented in 
SPiM and examples of models with varying structures, we refer to ll22l l5l[3l[ni4l[T9]|. 

3 Causality as Partial Order 

The stochastic semantics of the n calculus models in SPiM are given with continuous time Markov 
chains (CTMC). These CTMCs are infinite sets that collect the possible simulation trajectories of the 
modelled systems. In each possible trajectory of a process model, the duration of each state change is 
delivered by an exponential distribution from a continuous interval. Then, a simulation with a model is 
essentially a sequence of state transitions of the underlying model, where each transition is stamped with 
a time constant that imposes a total order on these transitions. In order to see this on an example, let us 
first define formally the class of models that we discuss in this paper. 

Definition 4 (model) Biochemical species are denoted with A,B,P, Q,. . .. A model is a pair (J^,^) 
where J" is a multise^of species denoting the initial state of the model and ^ is a finite set of chemical 
reactions of the form 

r:A^pPi + ...+Pi, or r :A + B Pi + ...+Pt 
such that k £N, r is the name of the rule, and p G is the rate of the reaction. 

Example 5 Consider the model with the reactions below and the initial state = ^A^, which can be 
implemented in SPiM as depicted in Figure |2] On the left-hand-side in Figure |j] we have a simulation 
output with this model. 

ri: A ^1.0 P + P r2:P^i.oB r^-.P^i.oC H-B + C^i.qD 

^Multisets are denoted by the curly brackets "\ |}". U , — and C denote the multiset operations corresponding to the usual 
set operations U , — and C , respectively. 
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directive sample 1000.0 new ch@r4 : chan ( ) and B() = ?ch; D() 

directive plot A() ; P(); and C() = !ch; () 

B(); CO; DO let A() = delayOrl; (P() I P()) 

and DO = O 

val rl = 1.0 val r2 = 1.0 and P() = do delay@r2; BO 

val r3 = 1.0 val r4 = 1.0 or delay@r3; CO run 1 of AO 



Figure 2: A SPiM implementation of the model given in Example |5] 

A simulation on such a process model can be seen as the interleavings of the computations of a 
place/transition Petri net, which can be infinite, e.g., in the case of polymerisation models HItI: in Petri 
nets, each transition has incoming arrows from places, containing tokens that are consumed as a result 
of the transition, and it has outgoing arrows to places, containing the tokens that are produced as a result 
of the transition. Then, each 'delay' action is a transition with a single incoming place, and each action 
resulting from the interaction of input ('?') and output (' ! ') actions is a transition with two incoming 
places for the input and output actions. For the purpose of this paper, we restrict the models to those, 
which can be expressed as finite sets of chemical reactions. This restriction allows us to denote each 
model interchangeably as a finite set of chemical reactions as well as the corresponding Petri net or the 
(equivalence class of the) SPiM implementation (since there does not exist a unique implementation for 
each set of reactions in SPiM). 

Example 6 Consider the Petri net depicted in Figure which is the net underlying the model given in 
Example^and implemented in Figure^ 

The consideration of a simulation trajectory of a process model from the point of view of Petri nets 
remains in agreement also with the CTMC semantics that imposes a total order on the transitions (see, 
e.g., 1201 ). However, when these transitions are inspected from the point of view of their dependencies 
on one another with respect to their production/consumption relationships, it is possible to relax this total 
order into a partial order by abstracting away from the time stamps of the state transitions. 

Example 7 Consider the model given in Example ^ In a simulation on this model, the reactions r2 
and rs are totally ordered due to the time stamps on them. However, we can relax this total order into a 
partial order with respect to the production/consumption relationship between the transitions as depicted 
in Figure^ 

In fact, the underlying Petri net of a process model hints this partial order structure on the transitions, 
as it can be seen on the example net depicted in Figure [4] However, the possible conflicts in simulations 
between various transitions with respect to the availability of resources remain implicit in the structure 
of a net. For example, consider a simulation state with the Petri net in Figure [4] where there is a token 
in place P, and thus there is a conflict between the transitions r2 and r^,. Although this conflict can be 
observed by inspecting multiple outgoing arrows from places, it is not always explicit in the structure of 
more complex Petri nets. 

Partial orders reflecting dependencies in computations have been studied as non-interleaving models 
of concurrency: in this setting, the synchronisation of concurrent transitions is realized by their com- 
mon successors and predecessors, which is in agreement with the notion of state transitions in Petri 
nets. When there are no resource conflicts between two partially ordered threads, they can take place 
in parallel (but not necessarily simultaneously), and their common predecessors and successors provide 
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a synchronisation mechanism. This is because the common predecessors deliver the resources for the 
parallel threads, and the completion of the execution of the threads delivers the resources required by 
their common successors. This view also provides a model of causal flow in the execution. 

A certain model of concurrency, which is closely related to Petri nets ll2ll . namely event structures, 
captures these ideas: in an event structure, each instance of a state transition is an event. The dependen- 
cies of events of a concurrent system are given with a partial order, and non-determinism is captured by 
a conflict relation, which is a symmetric irreflexive relation of events. When two events are in conflict, 
the execution of one of them instead of the other determines a different state space ahead. 

Definition 8 (labelled event structures) A labelled event structure (LES) is a structure {E,<,#,J^',£), 
where 

(i.) E is a set of events; 

( ii-) < ^ is a partial order such that for every e € E the set {e' £ E\e' < e} is finite, • 

(Hi.) the conflict relation # C is a symmetric and irreflexive relation such that ife#e' and e' < e" 
then e#e" for every e, e', e" G E; 

( iv.) ^ is a set of labels; 

(v.) £ : E ^ ^ is a labelling function. 

There are standard techniques in the literature for obtaining a labelled event structure from transition 
systems and from other models of concurrency ll27l . In particular, in flSl we present an application of 
these techniques for obtaining a labelled event structure from the multiplicative exponential linear logic 
encodings of Petri nets. We can analogously carry these ideas to the setting of the process models: for 
every model we use the operational semantics of the stochastic 7i calculus, given with the reduction 
rules in Definition|2j as a transition system. This way, we apply the techniques presented in lUSl to obtain 
the LES[[^]1 for the model 

Example 9 We obtain the event structure depicted in Figure^^or the model in Example^ 

Definition 10 Given a LES {E,<,#,^,i), for any event e G [ej denotes the set {e' € E\e' < e} of 
causes of event e. We say that 'i^ C E is a configuration if and only if (i.) for all e £'io we have that 
[eJ C '^; (ii.)for all e, e' G it is not the case that e#e'. 

Example 11 Consider the event structure given in Figure^ The sets {ei , 62, 63, 64} and {ei , 62, ej} are 
conflgurations. 



0. — > AO 

0.460566394485 A() — > P() P() 
0.734200976191 P() — > C() 
1.13404290455 P() — > B() 
1.86405154022 B() C() — > D() 




Figure 3: Graphical representation of the extraction of a partial order that reflects the produc- 
tion/consumption relationship of the transitions in a simulation with the model in Example [7] 
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Figure 4: The Petri net of the model given in Example [5] and the LES of this model, where the subscripts 
of the events coiTespond to the subscripts of the reactions(transitions) that they refer to. The dashed lines 
denote the conflict relation #. 

In [TSl, we have given an algorithm for extracting configurations from transition histories of Petri 
nets. Below, we analogously carry these ideas to the setting of the process models being discussed here: 
the CTMCs in SPiM treat the processes as resources in the sense that processes representing species are 
consumed and new ones are produced as a result of a reduction with SPiM. 'Resource consciousness' 
is the crucial ingredient that delivers the LES semantics in [181. We can thus extract a partial order that 
delivers the causal dependencies of the transitions of a particular simulation as sketched in Figure [3] Let 
us first collect some definitions that we need to formally describe this procedure on SPiM simulations. 

Definition 12 (simulation state) Let ^ = ( J^, ^) he a model with the reactions ri , . . . , r,j. A simulation 
state of denoted with 2^, is a multiset of species such that each species in is parameterised with a 
label k G denoting an id, a label 7?,^ G {init, ri, . . . , r„} denoting the reaction that created this species, 
and t G denoting its creation time. The initial simulation state 3fofor the initial state ^ of a model 
is obtained from J" by assigning a unique label k G to each species and marking each species with 
in it and the time 0. 

In SPiM, we can implement the notion of a simulation state by decorating each process with an 
integer valued id, which is communicated during the reduction of the process to its continuation. 

Example 13 Consider the SPiM program given in Figure |5] which is obtained from the program in 
Figure^by assigning an id to each process. Here, we have initially four instances of the process A. 
Thus, = H A(l, init,0), A(2, init,0), A(3, init,0), A(4, init,0) H is the initial simulation state for the 
initial state of the model implemented in Figure^^ 

Definition 14 (simulation trajectory) Let ^ = ( J^, J7) be a model. A simulation trajectory y on ^ 
is a list of structures of the form (t, L, R), called transition instances, where t G M+ is the time stamp of 
the transition. L and R are multisets of species that we call consumed and produced species, respectively. 
We assume that for each pair (L, R) there is a reaction r in 3" and the species on the left-hand-side and 
right-hand-side of x, respectively, are exactly those in L and R. In this case, we say that (t, L, R) is an 
instance of r. Vfe use [] to denote the empty trajectory and the operator '• ' to denote the concatination 
of a transition with a simulation trajectory, e.g., = (t, L, R) • 5^'. 

We use the simulation trajectories for representing the outcome of simulations with SPiM models. 
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directive sample 1000.0 
directive plot A() ; P(); 

BO; CO; DO 



val rl 
val r3 



1.0 
1.0 



val r2 
val r4 



1.0 
1.0 



let A(id:int) = delayOrl; 

(P(id) I P(id)) 

and P(id:int) = do delay@r2; B(id) 
or delayers ; C(id) 

and B(id:int) = ?ch; D(id) 



and C(id:int) = !ch; () 

and D(id:int) = () 

run ( A(l) I A(2) I 
A(3) I A(4) ) 



new ch(3r4 : chan O 



Figure 5: A SPiM implementation of the model given in Example [5] where each process is assigned an 
id as a parameter that is communicated during reductions. 



0. — > A(4) A(3) A(2) A(l) 
0.362742117504 A(4) — > P(4) P(4) 
0.403041278461 P(4) — > B(4) 
0.533386757223 P(4) — > C(4) 
1.38225471287 A(l) — > P(l) P(l) 
1.45347041154 A(2) — > P(2) P(2) 
1.51272758278 A(3) — > P(3) P(3) 
1.52969767363 P(3) — > B(3) 



1.53959162309 P(2) — > B(2) 

1.61357605349 P(3) — > B(3) 

1.62828014739 P(l) — > B(l) 

1.74376163407 P(2) — > C(2) 

1.93175167917 P(l) — > C(l) 

2.09975124488 B(4) C(l) — > D(4) 

2.10580589848 B(3) C(4) — > D(3) 

2.14571167381 B(2) C(2) — > D(2) 



Figure 6: A SPiM output of a simulation with the model given with the code in Figure |5] 

Example 15 Consider the simulation output given in Figure We denote this output as the simulation 
trajectory .5^ , where 5^' is a simulation trajectory that contains the transitions after the time point 1 .0. 

y= (0.362JA(4)^JP(4),P(4)^) . (0.403 J P(4) ^ J B(4) ^) . (0.533 J P(4) ^ J C(4) ^) . ^' 

Definition 16 (simulation configuration) Let ^ = (^, ^) be a model with the reactions ri, . . . , r<.. 
The simulation configuration/or a simulation trajectory 5^ and an initial simulation state 3^q, denoted by 
^y,%> the structure IJ-{^, 3^q), where the recursive junction }X on the pairs of simulation trajectories 
and simulation states is defined as follows. 

• if.y = (t, L, R) • .y such that L = \A{m) ^, R = -fl Pi (m), . . . ,Pk{m) ^, 

— there is a reaction r G {ri, . . . , r^;} such that (t, L, R) is an instance ofr, and 

— for transition label p and time label t', A(m,p,t') £ 3f, 

then }x{y,^) = {((m, p,t'),{m, r, t))} 

U At(^',(ir-|A(m,;.,0&)U|Pi(m,r,t),...,P,(m,r,t)^). 

• ify = {t,l,R)»y" such that L = {| A(mi),B(m2) ^, R = ^Pi{x), . . . ,Pt{x) \^,x = mi Vx = m2, 

— there is a reaction r G {ri, . . . , r^;} such that (t, L, R) is an instance ofr, and 

— for transition labels p,q and time labels t',t", A{mi,p,t') G 3f and A{m2,q,t") G 2f, 
then ^) = {{{mu p,t'),{x, r,t)), {{m2, q,t"),{x, r,t))} 

U At(^', {3f-Uim,P,t'),Aim2,q,t") ^) U I Pi(x, r,t), . . . ,P,(x, r,t) ^). 



Otherwise n{y,3f)=(d. 



28 



Flux Analysis in Process Models via Causality 



-• init 



0.362 




Figure 7: The simulation configuration obtained from the simulation trajectory delivered by the simu- 
lation output in Figure [6j and the graphical representation of the flux configuration obtained from it by 
merging the events that are instances of the same reactions. On the left, each node (x, r, ,f) is abbreviated 
with r; On the right, the thickness of the arrows are proportional with the weights of the edges. 

Example 17 Consider the model ^ with the SPiM code in Figure |^ the initial simulation state S^q 



given in Example 13 and the simulation trajectory 5^ given in Example 15 We get the simulation 



configuration Cy\% depicted in Figure^ 

The following proposition is a special case of a more general discussion presented in ifTSll . 



Proposition 18 Let ^ he a model with an initial simulation state 3fQ and a simulation trajectory 5^. 
The transitive reflexive closure of the simulation configuration C s^^sh <^ configuration in LES| 



4 Quantitative Flux Analysis on Models 

A simulation configuration contains the information on the causal dependencies of the transitions of a 
simulation. These causal dependencies are expressed as a partial order of individual events. However, 
despite the partial order representation of the causal flow, it remains possible to read the total order of 
the transitions given in the simulation with respect to the time stamps of these events. 

In this respect, a simulation configuration is a representation of a simulation, where the previously 
implicit information is made explicit. Although such a representation is rich in information, it is diffi- 
cult to make observations about the quantitative behaviour of the modelled system. Thus, a simulation 
configuration can be seen as raw data that needs to be processed further for analysis. 
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The idea here is to apply various transformations to the partially ordered structures of simulation 
configurations to extract useful information on the dynamics of the system. There are a variety of possi- 
ble choices for such transformations. In the following, we consider a transformation, which can provide 
a means for quantitative flux analysis in the stochastic setting of process algebra models. In this trans- 
formation, we merge the events that are instances of the same reaction to the same node, while counting 
the number of edges that collapse to a single edge. This way, we assign weights to the edges of the graph 
and obtain a cyclic graph with weights. 

Definition 19 (flux configuration) Let C y,% be a simulation configuration. The flux conflguration of 
C^,jj, denoted by fy.,^^, is the set that contains the triples of the form {p,q,n) such that {p,q,n) € 
F y^^g if and only if there are n number of pairs ((/, p, t),(j, q, t')) G C y,%for all the possible i,j,t,t'. 

Example 20 Consider the simulation configuration C depicted on the left-hand-side in Figure 
We obtain the flux configuration fy\% = {(init, ri,4), (ri, r2,5), (ri, r3,3), (r2, r4,3), (rs, r4,3)} which is 
depicted as the graph on the right-hand-side in Figure^ where the numbers on the edges are the weights 
of the edges. 

In comparison to the simulation trajectories and the simulation configurations, the flux configurations 
provide a more intelligible representation of the causal flow in a simulation from a quantitative point of 
view. In particular, in a flux configuration, the quantitative information on the flux of the species between 
reactions in the system is highlighted with respect to the simulation. This way it becomes possible to 
make observations on the system that are impossible, for example, by counting the number of species or 
reactions. However, it is important to note that there can be various choices to be made while processing a 
simulation configuration. For example, when we consider the transformation proposed above, we could 
have considered only causally independent events while merging them. While this criteria would not 
alter the graph in Figure [TJ it would deliver different graphs for the model that we discuss in the next 
section. 

5 A Biological Model: Rho GTP-binding Proteins 

In this section, we apply the ideas discussed above to a model of a biological system, that is, Rho GTP- 
binding proteins II17I [141 [5j| . 

The Rho GTP-binding proteins constitute a distinct family within the super-family of Ras-related 
small GTPases with twenty-two identified mammalian members, including Rho, Rac and Cdc42 lITTll . 
Rho GTP-binding proteins serve as molecular switches in various subcellular activities, regulating a va- 
riety of cell functions, including actin organisation and cell shape, cell adhesion, cell motility, membrane 
trafficking and gene expression 121. Their role can be perceived as regulating the transmission of an 
incoming signal further to effectors in a molecular module by cycling between inactive and active states, 
depending on being GDP or GTP bound, respectively. As depicted in Figure [8j GDP/GTP cycling is 
regulated by guanine nucleotide exchange factors (GEFs) that promote the GDP dissociation and GTP 
binding, whereas GTPase-activating proteins (GAPs) have the opposite effect and stimulate the hydrol- 
ysis of Rho-GTP into Rho-GDP. In the active GTP-bound state, Rho proteins interact with and activate 
downstream effectors. 

In (14], Goryachev and Pokhilko give an ordinary differential equations (ODE) model of the Rho 
GTP-binding proteins. The structure of this ODE model is depicted on the right-hand-side in Figure 
[8j the three forms of the Rho protein (GDP-bound RD, GTP-bound RT, and nucleotide free R) in the 
middle layer form complexes with GEE (E) in the bottom layer and with GAP (A) in the top layer. All 
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Plasma membrane 




Figure 8: Left: Rho GTP-binding protein cycle, adapted with permission from Macmillan Publishers 
Ltd: Nature 111], copyright 2002. Right: the structure of the Rho GTP-binding proteins model given in 
lfT4l . In the graph, R denotes the Rho GTP-binding protein, whereas RD and RT denote its GDP and 
GTP bound forms respectively. A and E denote GAP and GEF, respectively. Thus, RDE, for example, 
denotes the protein complex formed by RD and E. 



the reactions except GTP hydrolysis (2 — )• 1, 6 — 8, 5 — )• 3) are reversible. In their model, Goryachev and 
Pokhilko use the quantitative biochemical data on the Cdc42p member of the Rho family proteins. This 
results in an explanation of the experimentally observed rapid cycling of the Rho GTP-binding proteins 
between their GDP-bound off states and GTP-bound on states while displaying high activity, measured 
by the relative concentration of the GTP-bound Rho proteins (RT in Figure [8]l. The authors argue that 
GEF and GAP play distinct and separable roles in cycling control: the activity of Rho proteins is mainly 
delivered by the activity of the GEF and the turnover rate is a function of the GAP concentration. There- 
fore, to achieve a high activity and turnover rate simultaneously, the concentrations of GEF and GAP 
should be tightly controlled. Moreover, at large Eq and Aq, only a subset of the fluxes]^ of the original 
model is significant, while the remaining fluxes have negligible values. To test this hypothesis, Gory- 
achev and Pokhilko introduce a reduced model and provide a comparison with the original model with 
respect to the flux vectors that substantiates the claim. Goryachev and Pokhilko argue that in the efficient 
regime the operation of the GEF-GAP control module is based on two linear pathways: the GEF arm 
I— ^3— >-4— s-S— ;'2— s-l and the GAP arm 2 — 6 — 8 — 1, which are indicated by solid arrows in 
Figure [8] 

In [5|, Cardelli et al give a stochastic n calculus model of the Rho GTP-binding proteins, which 
is based on the ordinary differential equations model of (14). The model in |i51, displays an excellent 
agreement with the ODE model of [14] with respect to the RT activity on simulations with varying 
model structures and different regimes of initial concentrations. In the following, we consider this model 
for flux analysis by employing our software tool that implements the ideas above. 



■'in [141, the flux Ji,„ for a reaction that connects species / and m is defined with respect to the species concentrations and the 
rate constants, e.g., flux connecting RD and RDE is ^13 = A:i3.RD.E — ^31. RDE where ^13 and ^31 are the corresponding 
reaction rates. 
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Figure 9: A simulation witli tlie reduced model, given with the mid layer of the model depicted in Figure 
[8] and the graph representation of the flux configuration of a simulation for the interval <t < 140. 

5.1 Rho without Regulators 

At a first step for flux analysis with respect to simulations on a stochastic n calculus model of the Rho 
GTP-binding proteins, we consider the mid-layer of the model depicted on the right-hand-side of Figure 
[8] The SPiM code of the complete model is given in Appendix A. For the simulations on the reduced 
model, we run simulations on the complete model where the initial number of species for A and E are set 
as 0. The initial number of R is set as 1000. Because a variation on R levels does not cause a significant 
variation in the simulation behaviour as it can be observed from the results in (14] and |T|, here we do 
not consider variations on the Rq. The plot that displays the evolution of the number of species at a 
simulation that is run until the time t = 140 min. is depicted|^on the left-hand-side of Figure 9] 

The graph representation of the flux configuration for the time interval <t < 140 of a simulation is 
given on the right- hand-side of the Figure |9] In this plot, we observe that after the start of the simulation, 
the reactions R — RD and R — RT are competing for R and the flux ratio of these reactions reflect the bias 
of the rate constants of these reactions towards the reaction R — )• RT. Because the flux configuration is for 
the simulation until the time t = 140, the flux of the resources demonstrates the system's accumulation of 
steady state numbers of species from the start of the simulation: the initial concentrations are Rq = 1000, 
RDo = and RTq = 0. At the end of the simulation we have R = 1, RD = 503 and RT = 496. When we 
examine each reaction in the graph with respect to the weights of their incoming and outgoing arrows, we 
observe that the quantities correspond to those accumulated at the end of the simulation as demonstrated 



*The time unit of the model in 1141 is given in minutes. 
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Figure 10: Graph representation of the flux configuration of the simulation for the interval 100 <t < 140. 
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Figure 11: Graph representations of the flux configurations of the simulations for the intervals 100 < 
t < 120 and 120 < r < 140. 

below and they sum up to 1000. 

RT^RD: 1583- 1118 = 465RD ] RD ^ R : (1118 + 93) - (35 + 1175) = 1 R 1 

> sn^RD > 1 R 

R^ RD: (35 + 44 + 52)-93 = 38RD j RT ^ R : 1612 - (1560 + 52) = OR J 

R^ RT: (1560 + 956 + 1175) -(1583 + 1612) =496 RT 

By relying on the data that the cycling period of the system is 100 min [14], we analyse the flux 
configuration on the same simulation for the time interval 100 < f < 140. This provides an examination 
of the steady-state behaviour of the system. The graphical representation of the flux configuration for this 



interval is given on the left-hand-side of Figure 10 On the right-hand-side of this Figure, we simplify this 
graph by taking the sum of the fluxes for the cases, where there are opposite directional fluxes between 
two reactions. For example, we reduce the two arrows between the reactions R — ?■ RT and RT — >^ R to a 
single arrow. Below we examine the incoming and outgoing fluxes for each reaction. 

RT->RD:379-389 = -10RD RD ^ R : 395 - (389 + 6) = OR r ^ rj • 395 - (379+ 14) = 2 RT 
R^ RD: 14-6 = 8RD RT ^ R : 14- 14 = OR 

The observation that the fluxes of the individual reactions do not add up to reflects the fluctuations 
in the stochastic simulation. However, it is important to observe the sum — 10 + 8 + + + 2 = due 
to the fact that the simulation is performed on a closed system. In Figure [TT] we examine the flux 
configuration of the interval 100 < f < 140 further by separating it into two intervals 100 < f < 120 and 
120 <t < 140. We observe that in these smaller intervals, the fluxes of the reactions are closer to 0. 

5.2 Rho with GEF and GAP 

The behaviour of the model with the regulators GEF (E) and GAP (A) in lfT4l is insensitive to the R levels 
in terms of activity given by the RT/Rq ratio at the steady states of the simulations. The model has an 
activity maximum with the initial concentrations R = l.OmM, E = 0.776/mM and A = 0.66/xM. 

In order to analyse the flux of the system at the high activity regime, we run simulations with Rq = 
1000, Eo = 776 and Aq = 1, where we take the closest positive integer number value for Aq so that 
factoring of the other simulation parameters would not be required. This results in a near-maximum 
activity of approx 0.8 at the steady state with fluctuations due to stochastic simulation. A simulation plot 



with the model in Appendix A with these parameters is given on the left-hand-side of Figure 12 

We analyse the steady state behaviour of the system with respect to the simulation at this regime. 
The output of our tool, which displays the flux configuration of the simulation for the time interval 



KahramanoguUari 



33 




2.5 1-8 0,4 



Figure 12: Simulation plots with the complete model where the initial numbers of the species are Rq = 
1000 and Eq = 776. From left to right, the Aq value is set as 1, 10 and 100. An increase in Aq in the 
simulations causes decrease in the RT activity while reducing the recovery time. 



2.0 < t < 2.5, is given in Appendix B. From this flux configuration, we obtain the top-most graph in 
Figure [13] where we omit the fluxes which are less than 10% of the mean of the fluxes. This delivers the 
more dominant flux pathway in the system at steady state, given with 2— ;'6— s-S— s-l— s-S— ;'4— s-S— )-2, 
which is in agreement with the results of [14] with the exception that the path 2 — )• 1 is not included in 
our graph. The presence of this flux in the pathway could be explained only with the presence of a flux 
from the reaction RTE — RT to the reaction RT — RD. Although this flux can be read with a weight 
of 9 in our flux configuration ("22 ==> 9 appears 9 times ." in Appendix B), it is much smaller in 
weight in comparison to other fluxes that are depicted in the top-most graph in Figure 13 Moreover, we 
observe that the reaction RTE — RT + E and its reverse reaction RT + E — )• RTE deliver a strong cyclic 



flux, which cancels itself, thus they are not included in the graph in Figure 13 



When we carry this analysis to the simulations where Aq is increased to 10 and 100, we get the sim- 
ulation plots depicted in the middle and right-hand-side of Figure 12 We observe the same flux pathway 
patterns for these simulations as depicted as the middle and bottom graphs in Figure 13 However, we 
observe that an increase in Aq does not only decrease the RT activity, but also increases the turnover rate. 

It is important to observe that the graphs depicted in Figure 13 are generated from the model, and 
they reflect the predictions of the results in [141. Along these lines, Goryachev and Pokhilko argue 
that a subset of the fluxes of the original model is significant in the actual biological system, while 
the remaining fluxes have negligible values. Goryachev and Pokhilko substantiate this prediction by 



comparing the reduced model with the original model. In this respect, the graphs given in Figure 13 



depict a reduced model and agree with the predictions of II141I (with the exception of 2 — )• 1 as explained 
above). However, the graphs in Figure [13] result from the application of the ideas presented in Section 
[3] and Section |4] to the simulations, thus they are generated automatically, and their deUvery does not 
require an analysis of the model or the biological system being modelled. 



6 Discussion 

We have shown that the event structures view of simulations discussed above provide the means for flux 
analysis for process algebra models in a biologically meaningful manner. However, the ideas presented 
here should provide a starting point for further exploration rather than being conclusive. 

The relationship between process models and causality has been studied by various authors. In this 
regard, these previous investigations should provide touchstones for further investigations. Some of the 
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Figure 13: The graphs displaying the dominant fluxes with respect to the flux configurations obtained 
from simulations with Rq = 1000, Eq = 776 and varying initial Aq numbers. An increase in Aq increases 
the turnover rate, observed by comparing the ratio of the fluxes and the length of the time interval. 

works along these lines, which also provide further reference pointers, are as follows. In Danos 
et al draw connections between computational models of biological systems, event structures and their 
causality interpretation, while exploiting conflict as a mechanism of inhibition in a signalling pathway. 
The work by Pages and Soliman [12J attacks a similar problem by providing a formal interpretation of 
the relationship between reaction graphs and activation/inhibition graphs used by biologists. In lUl, Curti 
et al apply the ideas presented in iflOl to the n calculus models of biological systems where the causahty 
information on the modelled system is retrieved by labelling the syntax tree of the process expressions. 

Topics of future work include the application of the ideas above to BlenX [,23il language as well as 
an exploration of the relationship with the syntactical approach to causality with respect to enhanced 
operational semantics ifTOll . 
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Appendix A 

directive sample 4.0 

directive plot RD(); R() ; RT() ; 

RDEO; RE(); RTE() ; 

RDAO; RTAO; RA() 

val D = 50.0 
val T = 500.0 

new alQl . : chanO 
new a2@l .0:chan() 
new a3@l .0:chan() 

new e 1@0. 0054: Chan 
new e2@0. 0075: Chan 
new e3@0.43:chan() 



let AO = do ?al; () 

or ?a2; () 

or ?a3; () 

let E() = do ?el; () 

or ?e2; () 

or ?e3; () 



let R(id:int) = do 


delay(a0.033*D; RD(id) 


(* 


9 


-> 


1 


*) 


or 


delayOO.l*!; RT(id) 


(* 


9 


-> 


2 


*) 


or 


!a3; RA(id) 


(* 


9 


-> 


7 


*) 


or 


!e3; RE(id) 


(* 


9 


-> 


4 


*) 



and RD(id:int) = do 


delay©0.02; R(id) 


(* 


1 


-> 


9 *) 


or 


!al; RDA(id) 


(* 


1 


-> 


8 *) 


or 


!el; RDE(id) 


(* 


1 


-> 


3 *) 



and RT(id:int) = do 


delay@0.02; R(id) 


(* 


2 


-> 


9 


*) 


or 


delay@0.02; RD(id) 


(* 


2 


-> 


1 


*) 


or 


!e2; RTE(id) 


(* 


2 


-> 


5 


*) 


or 


!a2; RTA(id) 


(* 


2 


-> 


6 


*) 



and RE(id:int) = do 


delay@0.033*D; RDE(id) 


(* 


4 


-> 


3 


*) 


or 


delayOO.l*!; RTE(id) 


(* 


4 


-> 


5 


*) 


or 


delay(§1.074; ( R(id) I E() ) 


(* 


4 


-> 


9 


*) 



and RDE(id:int) = do delayOO. 136; ( RD(id) I E() ) (* 3 



-> 1 *) 
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or delayOe.O; RE (id) 




(* 


3 


-> 


4 


*) 


and RTE(id:iiit) 


= do delay(§76.8; ( RT(id) I 


E() ) 


(* 


5 


-> 


2 


*) 




or delay(§0.02; RDE(id) 




(* 


5 


-> 


3 


*) 




or delay(§0.02; RE (id) 




(* 


5 


-> 


4 


*) 


and RTA(id:int) 


= do delayOO . 0002 ; RA(id) 




(* 


6 


-> 


7 


*) 




or delay(§2 104.0; RDA(id) 




(* 


6 


-> 


8 


*) 




or delayOS.O; ( RT(id) I 


AO ) 


(* 


6 


-> 


2 


*) 


and RA(id:int) 


= do delay00.1*D; RDA(id) 




(* 


7 


-> 


8 


*) 




or delay(a0.0085*T; RTA(id) 




(* 


7 


-> 


6 


*) 




or delay(§500.0; ( R(id) I 


AO ) 


(* 


7 


-> 


9 


*) 


and RDA(id:int) 


= do delayOSOO.O; ( RD(id) 


1 AO ) 


(* 


8 


-> 


1 


*) 




or delayOO.02; RA(id) 




(* 


8 


-> 


7 


*) 



let ProduceR(id:int) = if id < 1001 

then (R(id) I ProduceR(id+l) ) 
elseO 



run 776 of E() 

run 1 of AO 

run 1 of ProduceR(l) 
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Appendix B 

The flux analysis output for the simulation with 
Ro = 1000, Eo = 776 and Aq = 1, for the time 
interval 2.0 <t < 2.5. 

Reactions : 



1 


A R 


— > RA fires times. 


2 


RA 


— > A R fires times. 


3 


RTA 


— > A RT fires times 


4 


R - 


-> RD fires times. 


5 


RD 


— > R fires times. 


6 


R - 


-> RT fires times. 


7 


RTE 


— > RDE fires 1 times. 


8 


RE 


— > RDE fires 2 times. 


9 


RT 


— > RD fires 9 times. 



10: RDE — > E RD fires 3 times. 

11: RE — > E R fires 2 times. 

12: E R — > RE fires 7 times. 

13: RTE — > RE fires 2 times. 

14: A RD — > RDA fires 16 times. 

15: A RT — > RTA fires 126 times. 

16: RTA — > RDA fires 126 times. 

17: E RD — > RDE fires 132 times. 

18: RDA — > A RD fires 142 times. 

19: RT — > R fires 6 times. 

20: RDE — > RE fires 149 times. 

21: RE — > RTE fires 145 times. 

22: RTE — > E RT fires 2240 times. 

23: E RT — > RTE fires 2095 times. 



Extracted dependencies : 



6 = 


==> 


23 


appears 


2 times. 


7 = 


==> 


20 


appears 


1 times. 


8 = 


==> 


20 


appears 


3 times. 


9 = 


==> 


14 


appears 


1 times. 


9 = 


==> 


17 


appears 


5 times. 


10 




17 


appears 


2 times . 


11 




12 


appears 


2 times. 


12 




21 


appears 


4 times. 


13 




21 


appeeirs 


1 times. 


14 




18 


appears 


16 times. 


15 




16 


appears 


126 times 


16 


==> 


18 


appears 


126 times 


17 




10 


appeeirs 


3 times. 


17 




20 


appears 


145 times 


18 




14 


appears 


15 times. 


18 




17 


appears 


125 times 


19 


==> 


12 


appears 


5 times. 


20 




8 


appears 


2 times. 


20 




11 


appeeirs 


2 times. 



20 


==> 


21 


appears 


140 times. 


21 




13 


appears 


1 times. 


21 




22 


appears 


148 times. 


22 




9 


appeairs 


9 times. 


22 




15 


appears 


126 times. 


22 




19 


appears 


6 times. 


22 




23 


appears 


2093 times 


23 




7 


appears 


L times. 


23 




13 


appears 


1 times. 


23 




22 


appears 


2092 times 



The flux analysis output for the simulation with 
Ro = 1000, Eo = 776 and Ao = 10, for the time 
interval 1.7 <? < 1.8. 

Reactions: 



1: 


RA — 


> RDA fires times. 


2: 


RTE - 


-> RDE fires times. 


3: 


RTA - 


-> A RT fires times. 


4: 


R — > 


RD fires times. 


5: 


RTE - 


-> RE fires times. 


6: 


RT — 


> RD fires times. 


7: 


RDE - 


-> E RD fires times. 


8: 


RD — 


> R fires 2 times. 


9: 


A R - 


-> RA fires 1 times. 


10 


RA - 


-> A R fires 1 times. 


11 


E R 


— > RE fires 10 times. 


12 


RE - 


-> RDE fires 3 times. 


13 


RT - 


-> R fires 1 times. 


14 


R — 


> RT fires 2 times. 


15 


RE - 


-> E R fires 10 times. 


16 


RTA 


— > RDA fires 101 times. 


17 


A RT 


— > RTA fires 101 times 


18 


A RD 


— > RDA fires 157 times 


19 


E RD 


— > RDE fires 123 times 


20 


RTE 


— > E RT fires 255 times 


21 


RDE 


— > RE fires 137 times. 


22 


E RT 


— > RTE fires 134 times 


23 


RDA 


— > A RD fires 261 times 


24 


RE - 


-> RTE fires 139 times. 



Extracted dependencies: 



7 ==> 18 appears 1 times. 

7 ==> 19 appeeirs 1 times. 

8 ==> 9 appears 1 times. 



KahramanoguUan 



8 = 


==> 


14 


appears 


times . 


9 = 


==> 


10 


appears : 


times . 


10 


==> 


11 


appeeirs 


1 times. 


11 


==> 


24 


appears 


13 times. 


12 


==> 


21 


appears 


5 times. 


13 


==> 


14 


appears 


1 times. 


15 


==> 


11 


appears 


9 times. 


16 


==> 


23 


appears 


102 times 


17 


==> 


16 


appears 


101 times 


18 


==> 


23 


appears 


159 times 


19 


==> 


21 


appears 


132 times 


20 


==> 


13 


appears 


1 times. 


20 


==> 


17 


appears 


101 times 


20 


==> 


22 


appears 


134 times 


21 


==> 


12 


appears 


3 times. 


21 


==> 


15 


appears 


10 times. 


21 


==> 


24 


appeeirs 


126 times 


22 




20 


appears 


117 times 


23 




8 


appears S 


times . 


23 


==> 


18 


appears 


156 times 


23 




19 


appeeirs 


122 times 


24 




20 


appeeirs 


138 times 



The flux analysis output for the simulation with 
Ro = 1000, Eo = 776 and Aq = 100, for the time 
interval 0.3 < ? < 0.4. 

Reactions : 



1 


RA 


— > RTA fires times. 


2 


RA 


— > RDA fires times. 


3 


R - 


-> RD fires times. 


4 


A R 


— > RA fires times. 


5 


RTA 


— > A RT fires times 


6 


R - 


-> RT fires times. 


7 


RT 


— > RD fires times. 


8 


RDA 


— > RA fires 1 times. 


9 


RE 


— > E R fires 1 times. 


10: RA 


— > A R fires 1 times. 


11: RD 


— > R fires 2 times. 



12 


E R 


— > RE fires 4 times. 


13 


RE - 


-> RDE fires 4 times. 


14 


RDE 


— > E RD fires 6 times. 


15 


E RT 


— > RTE fires 22 times. 


16 


E RD 


— > RDE fires 157 times. 


17 


RDE 


— > RE fires 143 times. 


18 


A RT 


— > RTA fires 148 times. 


19 


RTA 


— > RDA fires 149 times. 


20 


RE - 


-> RTE fires 140 times. 


21 


RTE 


— > E RT fires 168 times. 


22 


A RD 


— > RDA fires 2606 times 


23 


RDA 


— > A RD fires 2756 times 



Extracted dependencies : 



8 = 


==> 


10 


appears 


L times. 


9 = 


==> 


12 


appears 


L times. 


10 




12 


appeairs 


1 times. 


11 




12 


appears 


2 times . 


12 




20 


appears 


6 times. 


13 




17 


appears 


2 times. 


14 




22 


appears 


6 times. 


15 




21 


appeairs 


22 times. 


16 




14 


appears 


6 times. 


16 




17 


appears 


141 times. 


17 




9 


appears 


L times. 


17 




13 


appears 


4 times. 


17 


==> 


20 


appears 


134 times. 


18 




19 


appears 


149 times. 


19 




23 


appears 


151 times. 


20 




21 


appears 


146 times. 


21 




15 


appears 


22 times. 


21 


==> 


18 


appeeirs 


148 times. 


22 




8 


appears 


L times. 


22 




23 


appears 


2605 times 


23 




11 


appears 


2 times. 


23 




16 


appears 


157 times. 


23 


==> 


22 


appears 


2600 times 



